Caught while Dissolving: Revealing the Interfacial Solvation of the Mg2+ Ions on the MgO Surface

Interfaces between water and materials are ubiquitous and are crucial in materials sciences and in biology, where investigating the interaction of water with the surface under ambient conditions is key to shedding light on the main processes occurring at the interface. Magnesium oxide is a popular model system to study the metal oxide–water interface, where, for sufficient water loadings, theoretical models have suggested that reconstructed surfaces involving hydrated Mg2+ metal ions may be energetically favored. In this work, by combining experimental and theoretical surface-selective ambient pressure X-ray absorption spectroscopy with multivariate curve resolution and molecular dynamics, we evidence in real time the occurrence of Mg2+ solvation at the interphase between MgO and solvating media such as water and methanol (MeOH). Further, we show that the Mg2+ surface ions undergo a reversible solvation process, we prove the dissolution/redeposition of the Mg2+ ions belonging to the MgO surface, and we demonstrate the formation of octahedral [Mg(H2O)6]2+ and [Mg(MeOH)6]2+ intermediate solvated species. The unique surface, electronic, and structural sensitivity of the developed technique may be beneficial to access often elusive properties of low-Z metal ion intermediates involved in interfacial processes of chemical and biological interest.

1 Ambient Pressure Near Edge X-ray Absorption Fine Structure measurements  2 Processing of AP-NEXAFS spectra All measured Mg K-edge AP-NEXAFS spectra were subjected to background subtraction using the SNIP algorithm for background estimation, 1 and subsequently normalized dividing each background subtracted spectrum by its maximum intensity. Examples of raw experimental spectra, SNIP backgrounds and background corrected spectra are shown in Figure S2.
3 Decomposition of the AP-NEXAFS data into the spectra and fractional concentrations of key components Time-resolved spectroscopical measurements of chemical processes yield a series of spectra that may be positioned in a matrix D, where the columns of D are the spectra measured at time t. According to Lambert-Beer's law, at any given time a number N of "pure" and independent components weighed by their fractional concentration contributes to the measured signal. 2 Decomposing the experimental data into the spectra associated to the key species and in their relative concentration profiles can offer important insight in the investigated process. In the present work, such decomposition was performed with the PyFitit code, 2 a software that uses to such end an algorithm belonging to the MCR family. The starting point is the Singular Value Decomposition (SVD) equation: where the product U · Σ contains, on its N columns, a set of values that may be associated to the normalized absorption coefficients, Σ is a diagonal matrix known as the singular values term, whose elements are sorted in decreasing order, while V can be interpreted as the concentration matrix associated to the N-selected components. Lastly, the error matrix E represents the lack of fit between the experimental data matrix D and the reconstructed matrix µ = U ·Σ·V . The SVD deconvolution depends on the correct estimation of the number of components N present in the experimental spectral matrix. To this end, in this investigation we evaluated the percentage error committed in reproducing the experimental data with an increasing number N of components, as detailed in the main text, as shown in Figure S4. The percentage error function has been calculated with the following expression: where d ij and µ P C=n ij are the normalized absorbance values for the dataset and for the dataset reconstructed with N = n, respectively (K and m represent the number of acquired spectra and of the energy points, respectively, while n = 1, 2, ...K). At this point, all matrices in Equation 1 are solely mathematical solutions to the decomposition problem without physico-chemical meaning. Once N is established, the approach implemented by PyFitIt requires the introduction of a transformation N × N matrix T in Equation 1, using the relation I = T · T −1 : where the spectra belonging to the key species are given by S = U · Σ · T and their concentration profiles by C = T −1 · V . Subsequently, the matrix elements T ij of matrix T are modified by sliders to achieve S and C which are chemically and physically interpretable. Once this step is achieved, one can finally write: The unknown number of T ij elements of T is in principle equal to N 2 . In order to reduce such ambiguity, the AP-NEXAFS measured on the clean MgO surface was constrained to coincide with the first extracted spectral component for the analysis of the data related to the exposure of the MgO surface both to H 2 O and MeOH. This operation allows the reduction of the number of unknown T ij elements from N 2 to N 2 − N .
In the case of the H 2 O process, a 2 × 2 matrix T water was defined containing four elements. The solution for the decomposition presented in Equation 4 was obtained using the following matrix: Conversely, the following 2 × 2 matrix T methanol was defined, containing four elements, to decompose the UV-Vis data relative to the MeOH process:

Molecular dynamics simulations
Classical MD simulations were carried out on the Mg 2+ ion in aqueous solution and in methanol. Cubic boxes with one Mg 2+ ion and 500 water or MeOH molecules were built with side lengths of respectively 24.68 and 32.28Å, chosen in order to reproduce pure water and MeOH densities. Structures and interactions of the chemical species were represented with the SPC/E model for water, 3 while the OPLS force field was employed for MeOH. 4 Van der Waals interactions were introduced with the Lennard-Jones (LJ) potential with cross terms calculated with the Lorentz-Berthelot combining rules. The Mg 2+ ion was represented with LJ parameters optimized by Babu and Lim for the water case (σ = 2.430Å, = 0.0266 kcal mol −1 ). 5 A cut-off radius of 12Å was employed for all the non-bonded interactions and long-range electrostatic forces were taken into account with the Particle Mesh Ewald method. 6,7 Initial configurations were build with random positions with the PACKMOL package. 8 After energy minimization, each system was simulated in NVT conditions at 323 K for 10 ns, with the first 2 ns discarded as equilibration time. The temperature was controlled with the Nosé-Hoover thermostat with a coupling constant of 0.5 ps. Equations of motion were integrated with the leapfrog algorithm using a 1 fs time-step, while trajectories were saved every 100 steps. Stretching vibrations involving hydrogen atoms were constrained with the LINCS algorithm. 9 Simulations were performed with the Gromacs 2020.6 program. 10

DFT calculations
Minimum energy structures of clusters formed by the Mg 2+ ion with 4 and 6 water/MeOH molecules have been optimized at the density functional theory (DFT) level in gas phase. The B3LYP functional 11,12 was employed with the gaussian-type 6-31G(d,p) basis set for all the elements. Vibrational analysis was carried out to confirm the absence of imaginary frequencies and check that the stationary points were true minima. Calculations were carried out with the Gaussian 09 program. 13

Theoretical NEXAFS simulations
The Mg K-edge absoption spectra presented in this study were calculated using the FDMNES code, implementing the recently developed sparse solver method. [15][16][17] This software is based upon the Finite Difference Method (FDM), an attractive approach for the simulation of the photoelectron wave function beyond 100 eV above the absoption edge, avoiding the muffin tin approximation used in many common multiple scattering theory based codes. Specifically, the unit cell-normalized cross section σ(ω) was calculated as: where ω is the energy of the photon, α the fine structure constant, E g and E are the energies of the ground state Ψ (j) g and Ψ f , respectively, while the summation over j includes the contribution of all the atoms in the unit cells possessing index j. 18 The electron -photon interaction is treated classically employing the operator Θ, neglecting the magnetic part of the electromagnetic field and describing its electring portion with the first two terms of the multipolar expansion (corresponding to electric dipole and electric quadrupole excitations): where r is the relative position from the photoabsorber, is the photon polarization and k the photon wave vector. In all calculations the Schrōdinger-like equation was solved self-consistently to find the final states where there is a transition. 18 The calculated cross-sections were convoluted in a post-processing step by an energy-dependent arctangent function (Γ) in order to compare them to the experimental NEXAFS data. Γ is defined as follows: where E is the energy scale of the Mg K-edge NEXAFS spectrum, Γ i and Γ f are the core-level and final-state widths, respectively, E c and E w are the center and width of the arctangent function, respectively, while E f is the Fermi energy. 19 The density of electronic states of the investigated systems was also calculated through FDMNES. In order to reproduce the local environment of the Mg 2+ ion at the surface of MgO within a sufficient number of atomic planes, as probed by the AP-NEXAFS technique, the NEXAFS spectrum of MgO was simulated including in the calculation all scattering atoms within a cutoff radius of 5 A. The NEXAFS spectra obtained from the MD-extracted configurations were simulated including in the calculations water and methanol molecules within 6 A of the photoabsorber, i.e. until convergence was reached. A positive shift of 1294.9 eV was applied to all theoretical NEXAFS simulated spectra for the alignement with the experimental data.
7 Supplementary Figures S1-S9 Figure S1: a) 3D rendering of the reaction cell developed at APE-HE for the operando NEXAFS experiment. b) Side view of the reaction cell. Figure S2: Examples of the background subtraction procedure applied in this study to the raw measured Mg K-edge AP-NEXAFS spectra. The raw spectra, SNIP estimated backgrounds and background subtracted spectra are portrayed with blue, orange and green full lines, respectively.       The associated evolutions of the total absolute differences between averages of spectra computed with increasing N values are shown in the insets of panel.